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Abstract. - A thermodynamically consistent particle-based model for fluid dynamics with 
continuous velocities and a non-ideal equation of state is presented. Excluded volume inter- 
actions are modeled by means of biased stochastic multiparticle collisions which depend on 
the local velocities and densities. Momentum and energy are exactly conserved locally. The 
equation of state is derived and compared to independent measurements of the pressure. Re- 
sults for the kinematic shear viscosity and self-diffusion constants are presented. A caging and 
order/disorder transition is observed at high densities and large collision frequency. 



Introduction. - The efficient modeling of the long length- and time-scale dynamics of 
complex liquids such as colloidal and polymeric suspensions requires a simplified, coarse- 
grained description of the solvent degrees of freedom. A recently introduced particle-based 
simulation technique [1] — often called stochastic rotation dynamics (SRD) [2-7] or multi- 
particle collision dynamics [8, 9] — is a very promising algorithm for mesoscale simulations of 
this type. In additional to its numerical advantages, the algorithm enables simulations in the 
microcanonical ensemble, and fully incorporates both thermal fluctuations and hydrodynamic 
interactions. Furthermore, its simplicity has made it possible to obtain analytic expressions for 
the transport coefficients which are valid for both large and small mean free paths, something 
that is often very difficult to do for other mesoscale particle-based algorithms. This algorithm 
is particularly well suited for studying phenomena with Reynolds and Peclet numbers of order 
one, and it has been used to study the behavior of polymers [9, 10], colloids [1, 11-13], vesicles 
in shear flow [14], and complex fluids [15,16]. 

The original SRD algorithm models a fluid with an ideal gas equation of state. The fluid 
is therefore very compressible, and the speed of sound, c s , is low. In order to have negligible 
compressibility effects, as in real liquids, the Mach number has to be kept small, which means 
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that there are limits on the flow velocity in the simulation. It is therefore important to explore 
ways to extend the algorithm to model dense fluids. Our approach starts from what has been 
a common theme of most liquid theories, namely the separation of intermolecular forces into 
short- and long-range parts, which are then treated differently. The short-range component is 
a strong repulsion when molecules are close together; it leads to excluded volume effects which 
cause a decrease in the fluid's compressibility and eventual crystallization at low temperatures 
or high density. The long-range component is a weak attraction which can lead to a liquid- 
gas transition. The generic reference system for the short-range repulsive component of the 
force is the hard sphere system. In this letter we show how the SRD the algorithm can be 
modified to model excluded volume effects, allowing for a more realistic modeling of dense gases 
and liquids. This is done in a thcrmodynamically consistent way by introducing generalized 
excluded volume interactions between the fluid particles. The algorithm can be thought of as 
a coarse-grained multi-particle collision generalization of a hard sphere fluid, since, just as for 
hard spheres, there is no internal energy. In order to simplify the analysis of the equation of 
state and the transport coefficients, and enhance computational efficiency, the cell structure of 
the original SRD algorithm is retained. This work is a first step towards developing consistent 
particle-based algorithms for modeling, in the microcanonical ensemble, more general liquids 
with additional attractive interactions and a liquid-gas phase transition. 

Model. - As in the original SRD algorithm, the solvent is modeled by a large number N of 
point-like particles of mass m which move in continuous space with a continuous distribution 
of velocities. The system is coarse-grained into (L/a) d cells of a d-dimensional cubic lattice 
of linear dimension L and lattice constant a. The algorithm consists of individual streaming 
and collision steps. In the free-steaming step, the coordinates, r^t), of the solvent particles at 
time t are updated according to rj(t + r) = Ti(t)+TVi(t), where Vi(t) is the velocity of particle 
i at time t and r is the value of the discretized time step. In order to define the collision, we 
introduce a second grid with sides of length 2a which (in d = 2) groups four adjacent cells 
into one "supercell" . 

As discussed in Refs. [2] and [3], a random shift of the particle coordinates before the 
collision step is required to ensure Galilean invariancc. All particles are therefore shifted 
by the same random vector with components in the interval [—a, a] before the collision step 
(Because of the supercell structure, this is a larger interval than in the conventional SRD 
algorithm). Particles arc then shifted back by the same amount after the collision. To initiate 
a collision, pairs of cells in every supercell are randomly selected. As shown in Fig. 1, three 
different choices are possible: a) horizontal (with a\ = x), b) vertical {or 2 = y), and c) diagonal 
collisions (with cr 3 = (x + y)/V2 and er 4 = (x — y)/\/2). Note that diagonal collisions are 
essential to equilibrate the kinetic energies in the x— and y— directions. 

In every cell, we define the mean particle velocity, 



where the sum runs over all particles, M n , in the cell with index n. The projection of the 
difference of the mean velocities of the selected cell-pairs on aj, Au — aj ■ (ui — u 2 ), is then 
used to determine the probability of collision. If Au < 0, no collision will be performed. For 
positive Au, a collision will occur with an acceptance probability which depends on Au and 
the number of particles in the two cells, Mi and M 2 . This rule mimics a hard-sphere collision 
on a coarse-grained level: For Au > clouds of particles collide and exchange momenta. For 
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Fig. 1 Fig. 2 



Fig. 1 - Schematic of collision rules. Momentum is exchanged in four ways, a) horizontally along a\, 
b) vertically along (T2, c) diagonally and d) off-diagonally along cr 3 and (T4 respectively, according to 
Eq. @. w and Wd denote the probabilities of choosing collisions a), b) and c), d) respectively. 

Fig. 2 - Diffusion coefficient as a function of r. The data points in the inset are data for the shear 
viscosity measured using Green-Kubo relations, as a function of rksT. The solid line shows the 
analytical result, Eq. Parameters: L = 64a, M — 5, fcsT = f.O and A = f/60. 

reasons discussed in the following, we have used the acceptance probability 

p A (Mi,M 2 ,Au) = 8(Au) tanh(A) with A = iA«M 1 M 2 , (2) 

where O is the unit step function and A is a parameter which allows us to tune the equation 
of state. The hyperbolic tangent function was chosen in @ in order to obtain a probability 
which varies smoothly between and 1. 

Once it is decided to perform a collision, an explicit form for the momentum transfer 
between the two cells is needed. The collision should conserve the total momentum and 
kinetic energy of the cell-pairs participating in the collision, and in analogy to the hard-sphere 
liquid, the collision should primarily transfer the component of the momentum which is parallel 
to the connecting vector Uj. In the following, this component will be called the parallel or 
longitudinal momentum. There are many different rules which fullfill these conditions. Our 
goal here is to obtain a large speed of sound. We therefore use a collision rule which leads to 
the maximum transfer of the parallel component of the momentum and does not change the 
transverse momentum. The rule is quite simple; it exchanges the parallel component of the 
mean velocities of the two cells, which is equivalent to a "reflection" of the relative velocities, 

v!(t + r)-J=-(v!(t)-J), (3) 

where it" is the parallel component of the mean velocity of the particles of both cells. The 
perpendicular component remains unchanged, 

vt{t + r)=vi{t). (4) 

It is easy to verify that these rules conserve momentum and energy in the cell pairs. 
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Because of x — y symmetry, the probabilities for choosing cell pairs in the x— and in- 
directions (with unit vectors ct\ and <r 2 in Fig. 1) are equal, and will be denoted by w. The 
probability for choosing diagonal pairs (cr 3 and 0*4 in Fig. 1) is given by Wd — 1 — 2w. w and 
Wd must be chosen to that the hydrodynamic equations are isotropic and do not depend on the 
orientation of the underlying grid. This can be done by considering the temporal evolution 
of the lowest moments of the velocity distribution function. It is sufficient to consider the 
following three moments for a single particle i, 

( (vl(t)) 

\ (Vix{t)V iy (t)) 

Assuming molecular chaos, the collision rules can be used to determine the eigenvalues of the 
relaxation matrix, R, defined by *&i(t + r) = R^i(t). 

Because of the conservation of energy, one of the three eigenvalues of R is equal to one; 
the other two are given by Ai = Wd + 2w(2/M — 1) and A2 = 2w + Wd(2/M — 1), where M is 
the average number of particles per cell. Isotropy requires that Ai = A2, a condition that can 
be fullfiled for arbitrary M only if Wd = 1/2 and w = 1/4. Simulations show that both the 
speed of sound and the shear viscosity are isotropic for this choice. Note, however, that this 
does not guarantee that all properties of the model are isotropic. This becomes apparent at 
high densities or high collision frequency, 1/t» 1, where inhomogenuous states with cubic 
or rectangular order can be observed (see Fig. 4 and accompanying discussion). 

Transport coefficients. - The transport coefficients can be determined using the same 
Green-Kubo formalism as was used for the original SRD algorithm [2-4]. In particular, the 
kinematic shear viscosity is given by 

00 ' 

V= Nk T ^ ( S xy(0)S xy (nT)), (6) 

where 

N 

S X y{nr) = J2 {v jx (nT)A£ jv (nT) + Av jx {nr)[A^ y (nr) - z s jly ([n + l]r)/2]) (7) 
i=i 

is the off-diagonal element of the stress tensor S. and £?(t) are the cell coordinates of 

particle j in the fixed and shifted frames at time t, respectively, AfcJt) = + r) — £j(t), 
A£j (t) = £j (t + r) — £j (t + r), and Avj (t) = Vj (t + r) — Vj (t) . Zj t indexes pairs of cells which 
participate in a collision event; the second subscript, /, is the index of the collision vectors cri 
listed in Fig. 1. For example, for collisions characterized by cri, Zj lx — 1 if ^ x in Q is one 
of the two cells on the left of a supercell and z S j lx = — 1 if £J X is on the right hand side of a 
supercell; all other components of z s are zero. In general, the components of are either 0, 
1, or — 1. Using {z^/}, the collision invariants of the model can be written as 

3 

where aij = 1 for the density, {a^j} = {vp_ij} are components of the particle momentum, 
and ad+2,j = v j/2 is the kinetic energy of particle j [3]. The analogous collision invariants for 
the standard SRD algorithm are given in Eq. (25) of [3]. The vectors z 8 are constructed so 
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that the sum of the two exponentials in © is the same for two particles if and only if they 
are in partner cells in a collision with index I (see Fig. 1). 

The self-diffusion constant D is given by a sum over the velocity-autocorrelation function 
(see, e.g. Eq. (102) in [3]) and can be evaluated analytically assuming molecular chaos. Due to 
the excluded volume interactions, density fluctuations are supressed in the current algorithm; 
ignoring these fluctuations, one finds 



which is in good agreement with simulation data, see Fig. 2. 

Equation of state. - The collision rules conserve the kinetic energy, so that the internal 
energy of our system should be the same as that of an ideal gas. Thermodynamic consistency 
requires that the non-ideal contribution to the pressure is linear in T. As will be shown, this 
is possible if the coefficient A in J2J) is chosen small enough (see Fig. 3). 

We use here the mechanical definition of pressure — the average longitudinal momentum 
transfer across a fixed interface per unit time and unit surface area — to determine the equation 
of state. We consider only the momentum transfer due to collisions, since that coming from 
streaming constitutes the ideal part of the pressure. 

Take an interface that is parallel to the y— axis and consider the component p xx of the 
pressure tensor. Only collisions with label I = 1, 3, and 4 of the collision vector <T; in Fig. 1 
contribute to the momentum transfer in this case. Consider the contribution to the momentum 
transfer across the cell boundary from collisions with 1 = 1. For fixed number of particles, 
Mi and M 2 , in the two cells, the thermal average of the momentum transfer, AG X , across the 
dividing line is 



The factor 1/2 comes from the position average of the dividing line, since the collision occurs 
n the shifted cells, and the integral is restricted to positive Au because the acceptance rate is 
zero for Au < 0. pa(Au) is the probability that u\ x — U2 X for the micro-state of two cells is 
equal to Au. w — 1/4 is the probability of selecting this collision. 

Expanding the acceptance probability, Eq. (J2J) , in A = A Au Mi M 2 leads to 
Pa(Mi, M 2 , Au) = 0(Au)(A — A 3 /3 + ...). The contributions to the pressure from all terms of 
this series can be calculated, but since the resulting contribution to the pressure from a term 
proportional to A™ is of order T", we restrict ourselves to the first term. 

The resulting contribution to the pressure, P(tri, Mi, M 2 ), for fixed Mi and M2 is the 
average momentum transfer per unit area and unit time, so that using Eqs. I|10|) . we have 
P(cti,Ml,M 2 ) = wAk B T M 1 M 2 /(2aT) + 0(A 3 T 2 ). A similar calculation can be performed 
for the contributions from the diagonal collisions, which occur with the probability Wd- Using 
w = 1/4 and uid = 1/2 and averaging over the number of particles per cell (assuming that 
they are Poisson-distributed and that the particle number distributions in adjacent cells are 
not correlated), one finds the non-ideal part of the pressure, 



where Pid = k B T M/a 2 is the ideal gas contribution to the pressure (in d — 2). Note that the 
same result is obtained if, instead of averaging over Mi and M 2 , we simply set Mi = M 2 = M, 
the average number of particles per cell. P n is quadratic in the particle density, p = M/a 2 , 






(10) 




(11) 
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Fig. 3 
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Fig. 4 



Fig. 3 - Non-ideal pressure, P n , as a function of fcsT/r averaged over 10 5 time steps. Both ksT and 
r ranged from 0.005 to 4. The line represents the theoretical expression, Eq. dill) . For r = 0.005 and 
feT = 1, P n is five times larger than P i( i. Parameters: L — 64a, M = 5, A = 1/60. 

Fig. 4 - Freezing snapshot after 10 6 time steps. Parameters: L x — L y = 32a, M = 5, fcgT = 
3.125 x 10' 5 , r = 10~ 3 , A = 1/60. 



as one would expect from a virial expansion. The prefactor A must be chosen small enough 
that higher order terms in this expansion are negligible. We have found that prefactors A 
leading to acceptance rates of about 20% are sufficiently small to guarantee that the pressure 
is linear in T (see Fig. 3). In order to measure P n , we have used the fact that the average of 
the diagonal part of the microscopic stress tensor gives the virial expression for the pressure 

P = P u +P n = (j2 {'v.- + A« isB [AC jx - z s jlx /2] } ^ . (12) 

The first term, (v x .jA^ x j) = (tv% •), gives Pjd, as discussed in Ref. [3], The average over 
the second term vanishes (see Ref. [3]), while the average of the third term is the non-ideal 
part of the pressure, P n . Simulation results for P n obtained using $T2$ are in good agreement 
with the analytical expression, (|llf) (see Fig. |3J). In addition, measurements of the density 
fluctuations, {\pk\ 2 }, at small wave vectors k, as well as results for the adiabatic speed of 
sound obtained from simulations of the dynamic structure factor, are both in good agreement 
with the predictions following from Eq. (|llfl . These results provide strong evidence for the 
thermodynamic consistency of the model. 

Caging and order/disorder transition. If the non- ideal part of the pressure is large 
compared to the ideal pressure, ordering effects can be expected. For small A, both contribu- 
tions to the pressure are proportional to the temperature, so that just as in a real hard-sphere 
fluid, changing the temperature does not lead to an order/disorder transition. On the other 
hand, the two contributions to the pressure have different dependencies on the density and 
time step, r. In fact, r can be interpreted as a parameter describing the efficiency of the 
collisions; lowering r results in a higher collision frequency, and has a similar effect to making 
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the spheres larger in a real hard-sphere system. We therefore expect caging and ordering effect 
if either p is increased or r is decreased. This is indeed the case. For t < 0.0016, p — 5 and 
ksT — 1.0, an ordered cubic state is observed. The cubic symmetry of the ordered state is 
clearly an artifact of the cubic cell structure, and it would be interesting to see if this could be 
removed by using an hexagonal cell structure or incorporating random rotations of the grid. 
One of the surprising features of this crystalline-like state is, that x — y symmetry can be 
broken. Furthermore, there is the possibility of having several metastable crystalline states 
corresponding to slightly different lattice constants and number of particles per "cloud" . As 
expected, the lattice constants of these ordered states are slightly smaller than the super-cell 
spacing, 2a, which sets the range of the multi-particle interaction. In this state, the diffusion 
coefficient becomes very small; particles are caged and can barely leave their location. 

To understand this behavior, note that without collisions, particle clouds will broaden 
due to streaming; this will happen faster the higher the temperature. Due to the grid shift, 
particles at the perimeter of the clouds will more often undergo collisions with neighbor clouds. 
These collisions backscatter the particles, forcing them to fly back towards the center of their 
cloud. There is a correlation between the distance from cloud center and rate at which it is 
backscattered, leading to stable cloud formation. A particle which is left alone between clouds 
will feel repulsion from all clouds and moves around very quickly until it is absorbed into a 
cloud. 

Conclusion. - The model presented in this letter is the first extension of the SRD al- 
gorithm to model fluids with a non-ideal equation of state. It was shown that the model 
is thermodynamically consistent for the correct choice of acceptance probabilities and repro- 
duces the correct isotropic hydrodynamic equations at large length scales. Expressions for the 
equation of state and the self-diffusion constant were derived and shown to be in good agree- 
ment with numerical results. Simulation results for the kinematic viscosity were presented, 
and it was shown that there is an ordered state for large densities and collision frequencies. 
A detailed analysis of the transport coefficients will be presented elsewhere. 
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